% This code produces the main figure for the monetary part of the paper
% by Raphael schoenle, February 1, 2016 and April 2019
% edited by Jiacheng Feng modified March 18, 2019

% six cases with parameterizations

clearvars -except wd 

% load sector code labels
% initial data
 agg_1 = fopen('data_341_NAICS_sector_names.csv');
 fmt = '%s%s';
 data_1 = textscan(agg_1,fmt,'headerlines',1,'delimiter',',');
 agg_2 = fopen('gdpshares_58_2002_valad_concorded_new.csv');
 fmt = '%s%s';
 data_2 = textscan(agg_2,fmt,'headerlines',1,'delimiter',',');
% BEA NAICS codes: these are always 341 identical codes from the first column
 bea_codes_341 = data_1{1};
 bea_codes_58 = data_2{1};

% Data
load FPA.mat
%Parameters values
params.rho_a = 0.9;         % persistence agg tech shocks
params.rho_mu = 0.9;         % persistence monetary shock
params.rho_k = 0.9;         % persistence idiosync tech shock
params.beta = 0.9975;      % discount factor
params.theta = 6;          % elast of subst within sectors
params.delta = 0.5;

Ck = csvread("2002CkDetailed.csv");
Znum = csvread("2002RevShareDetailed.csv");
n = length(Ck); % Number of Industries
Z = Znum(1:n,1:n);

Wc = Ck/sum(Ck); % Share of consumption; the weights omega
W = Z./(ones(n,1)*sum(Z)); %Normalize the share of inputs for each sector, omega
psi = params.delta*(params.theta - 1)/params.theta;
psi_C = (1-params.delta*(params.theta - 1)/params.theta);
n_k = inv(eye(n)-psi.*W)*(psi_C*Wc);

%number of sectors
n_sectors = 341;
n_sectors1 = 58;

%number of periods in IRF
n_periods = 30;
n_periods_idio = 30;

%number of cases
num_case = 7;

%pre-definition of variables
% third dimension: number of shocks: monetary, technology, idiosyncratic
c_imp = zeros(n_periods, num_case, 3);
inf_imp = zeros(n_periods, num_case, 3);
mc_imp = zeros(n_periods, num_case, 3, n_sectors);
pos_akshock = zeros(n_sectors,1);
pos_akshock1 = zeros(n_sectors1,1);

% parameter sets: only base case for Figure 1!
delta_set = [0.5];
eta_set = [2];
phi_set = [2 ];
sigma_set = [1];
phi_pi_set = [1.24];
phi_c_set = [0.33]./12;
phi_rho_i_set = [0 ];
phi_gc_i_set = [0];
n_parameters = length(phi_rho_i_set);

% Position indices for shocks
pos_mushock = 1;
pos_ashock = 2;
for i = 1:n_sectors
    pos_akshock(i) = 2 + i;
end
for i = 1:n_sectors1
    pos_akshock1(i) = 2 + i;
end
     
% sectoral weights in consumption
Omega_C_het = Wc;
Omega_C_hom = 1/n_sectors*ones(size(Omega_C_het));

% 58 sectors
wc_58_het = xlsread('gdpshares_58_2002_valad_concorded_new.xlsx');
Omega_C_58_het = wc_58_het(:,2);
Omega_C_58_hom = 1/n_sectors1*ones(size(Omega_C_58_het));

% sectoral Calvo
% 350 sectors first
FPA_het  = FPA;
FPA_hom = mean(FPA_het)*ones(size(FPA_het)); 
FPA_flex = 0.99*ones(size(FPA_het));
FPA_sticky = 0.01*ones(size(FPA_het));
% Omega_C weighted mean
FPA_hom_weighted = Omega_C_het'*FPA_het*ones(size(FPA_het));

% 58 sectors
FPA_58_het = xlsread('secfreq_58.xlsx');
% adjust to have same mean as 350 sectors
FPA_58_het = (mean(FPA_het)-mean(FPA_58_het))+FPA_58_het;
FPA_58_hom = mean(FPA_het)*ones(size(FPA_58_het));
FPA_58_flex =  0.99*ones(size(FPA_58_het));
FPA_58_sticky =  0.01*ones(size(FPA_58_het));
FPA_58_hom_weighted = Omega_C_58_het'*FPA_58_het*ones(size(FPA_58_het));

%sectoral weights 
% Omega_hat; [rows = supplier, columns = customer]: input shares over all
% suppliers add up to 1
Omega_het = W;
Omega_hom = 1/n_sectors*ones(size(Omega_het)); 
Omega_C_as_IO = repmat(Omega_C_het',n_sectors,1)';

Omega_58_het = xlsread('IO_links_58.xlsx')';
Omega_58_het = (Omega_58_het./repmat(sum(Omega_58_het,2),1,n_sectors1))';
Omega_58_hom = 1/n_sectors1*ones(size(Omega_58_het)); 
Omega_58_C_as_IO = repmat(Omega_C_58_het',n_sectors1,1)';

%% 341 Sectors
%Case1: heterogeneous Calvo, heterogeneous size, heterogeneous I/O
FPA(:,:,1) = FPA_het;
Omega_C(:,:,1) = Omega_C_het;
Omega(:,:,1) = Omega_het;

%Case2: heterogeneous Calvo, heterogeneous size, homogeneous I/O
FPA(:,:,2) = FPA_het;
Omega_C(:,:,2) = Omega_C_het;
Omega(:,:,2) = Omega_hom;

%Case3: heterogeneous Calvo, homogeneous size, heterogeneous I/O
FPA(:,:,3) = FPA_het;
Omega_C(:,:,3) = Omega_C_hom;
Omega(:,:,3) = Omega_het;

%Case4: heterogeneous Calvo, homogeneous size, homogeneous I/O
FPA(:,:,4) = FPA_het;
Omega_C(:,:,4) = Omega_C_hom;
Omega(:,:,4) = Omega_hom;

%Case5: homogeneous weighted mean Calvo, homogeneous size, homogeneous I/O
FPA(:,:,5) = FPA_hom_weighted;
Omega_C(:,:,5) = Omega_C_hom;
Omega(:,:,5) = Omega_hom;

% for Appendix only
%Case6: heterogeneous Calvo, heterogeneous size, heterogeneous
%I/O based on size weights
FPA(:,:,6) = FPA_het;
Omega_C(:,:,6) = Omega_C_het;
Omega(:,:,6) = Omega_C_as_IO;

%Case7: flex-price Calvo, homogeneous size, homogeneous I/O
FPA(:,:,7) = FPA_flex;
Omega_C(:,:,7) = Omega_C_hom;
Omega(:,:,7) = Omega_hom;


% outer loop over the parameter values, index k
for k = 1:n_parameters
    
    params.eta = eta_set(k);
    params.frisch = phi_set(k);
    params.sigma = sigma_set(k);
    params.phi_pi = phi_pi_set(k);
    params.phi_c = phi_c_set(k);
    params.phi_gc = phi_gc_i_set(k);
    params.rho_i = phi_rho_i_set(k);
     
    % inner loop over the cases, index i
    % 1-6: 341 sectors
    % 7-12: 58 sectors
    
    % all sector cases
    for i = 1:7
    tic
    i
    
    % delete Table 5a,b output if existing
	if exist(['Table5a_case_' num2str(i-num_case/3) '_params' num2str(k) '_58.csv'], 'file')==2
        delete(['Table5a_case_' num2str(i-num_case/3) '_params' num2str(k) '_58.csv'])
    end
    if exist(['Table5b_case_' num2str(i-num_case/3) '_params' num2str(k) '_58.csv'], 'file')==2
    	delete(['Table5b_case_' num2str(i-num_case/3) '_params' num2str(k) '_58.csv'])
    end
    if exist(['Table5a_case_' num2str(i) '_params' num2str(k) '_341.csv'], 'file')==2
        delete(['Table5a_case_' num2str(i) '_params' num2str(k) '_341.csv']);
	end
    if exist(['Table5b_case_' num2str(i) '_params' num2str(k) '_341.csv'], 'file')==2
    	delete(['Table5b_case_' num2str(i) '_params' num2str(k) '_341.csv']);
    end
    
            % 341-sector cases
            %%Generate Gensys matrices
            [G0, G1, C, PSI, PI, pos_c, pos_pc, pos_mc, pos_yk, pos_p, pos_ps, pos_zk]  = ...
                Gensys_matrix(n_sectors, FPA(:,:,i),Omega_C(:,:,i), Omega(:,:,i), params);
        toc
            %%Call Gensys
            [G1, C, imp, fmat, fwt, ywt, gev, eu, loose] = gensys(G0, G1, C, PSI, PI);
        toc
            %check existence and uniqueness
            eu
        
            yt_imp_mu = Gensys_IRF(G1, n_sectors, n_periods, pos_mushock, imp);
            filename = ['temp' num2str(i) '.mat'];
            save(filename,'yt_imp_mu');
            
            % output = [index for sectors, consumption IRF, inflation IRF, real MC IRF]
            output = [1:n_periods; yt_imp_mu(pos_c,2:end); ...
                yt_imp_mu(pos_pc,2:end) - yt_imp_mu(pos_pc,1:end-1); ...
                yt_imp_mu(pos_mc,2:end) - (yt_imp_mu(pos_pc,2:end)) ]';
            autocorr_c = autocorr(output(:,end-1),1);
            autocorr_inf = autocorr(output(:,end),1);
            autocorr_mc = autocorr([mean(yt_imp_mu(pos_mc,2:end)) - (yt_imp_mu(pos_pc,2:end))]',1);
            % matrix size [n_periods x 7] ??
            
            % output matrix: initial C, Pi, MC, cumulative C, Pi, mean(MC)...
            % AR(1) of C, Pi, MC + IRFs (C, Pi, MC)
            output_matrix = [output(:,1) repmat([output(1,2:3)...
                mean(yt_imp_mu(pos_mc,2) - (yt_imp_mu(pos_pc,2)))...
            sum(output(:,2:3)) sum(yt_imp_mu(pos_mc,2:end) - (yt_imp_mu(pos_pc,2:end)))...    
            autocorr_c(2) autocorr_inf(2) autocorr_mc(2)],n_periods,1) output(:,2:end)];

            % cumulative responses - sort and get sector indices
            % use ck = yk - zk
            [output_top_bottom, I] = sort(sum(yt_imp_mu(pos_yk,2:end),2));
            top_bottom = [I(1:10) ; I(end-9:end)];
            %[output_top_bottom(1:10) output_top_bottom(end-9:end)]'

            % save cumulative response and identity of extreme bottom/top 10 sectors
            % output all sectors 
            filename = ['Table5a_case_' num2str(i) 'params' num2str(k) '_341.csv'];
            dlmwrite(filename, [output_top_bottom(1:10) output_top_bottom(end-9:end)]', '-append',  'precision', '%.10f', 'delimiter', '\t');

            filename = ['Table5b_case_' num2str(i) 'params' num2str(k) '_341.csv'];
            fid = fopen(filename, 'a');
            fprintf(fid, '"%s"\n', [ strjoin(bea_codes_341(I))]);
            fclose(fid);

            % save output
            filename = ['Gensys_mu_output_main_case_' num2str(i) 'params' num2str(k) '_341.csv'];
            fid = fopen(filename, 'w');
            fprintf(fid,'Period\tConsumption_impact_irf_mu\tInflation_impact_irf_mu\tMC_mean_impact_irf_mu\tConsumption_cum_irf_mu\tInflation_cum_irf_mu\tMC_mean_cum_irf_mu\tConsumption_persistence_mu\tInflation_persistence_mu\tMC_persistence_mu\tConsumption_mu\tInflation_mu\tMC_mu \n');
            fclose(fid);
            dlmwrite(filename, (output_matrix), '-append', 'precision', '%.10f', 'delimiter', '\t');
            filename = ['Table5c_case_' num2str(i) '_params' num2str(k) '_341.csv'];
            dlmwrite(filename, [I [1:length(I)]'] );
          
            % save output as matlab file for later plotting
            filename = ['Gensys_mu_output_main_case_' num2str(i) 'params' num2str(k) '_341.mat'];
            save(filename, 'output_matrix')

    
        % for plotting
        c_imp(:,i,1) = yt_imp_mu(pos_c,2:end)';
        inf_imp(:,i,1) = (yt_imp_mu(pos_pc,2:end) - yt_imp_mu(pos_pc,1:end-1))';
        mc_imp(:,i,1) = [mean(yt_imp_mu(pos_mc,2:end)) - yt_imp_mu(pos_pc,2:end)]';

    end
    
end

save('plot_data_fig1.mat','c_imp','inf_imp','mc_imp');
n_periods = 30;
n_sectors = 341;
k = 1;

    % make figures for each parameter combination k
    % 341 sectors
    i = 1
    figure('units','normalized','position',[.1 0 .4 1])
    % 341 goods, inflation
    subplot(3,1,1)
    plot(1:n_periods,c_imp(:,1,1)','-b','linewidth',2)
    hold on
    plot(1:n_periods,c_imp(:,i+1,1)','-b+','linewidth',2)
    plot(1:n_periods,c_imp(:,i+2,1)','--b','linewidth',2)
    plot(1:n_periods,c_imp(:,i+3,1)',':b','linewidth',2)
    plot(1:n_periods,c_imp(:,i+4,1)','-.b','linewidth',2)
    title('\textbf{Consumption}','Interpreter','latex')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
    % 341 goods, mean MC
    subplot(3,1,2)
    plot(1:n_periods,inf_imp(:,i,1)','-b','linewidth',2)
    hold on
    plot(1:n_periods,inf_imp(:,i+1,1)','-b+','linewidth',2)
    plot(1:n_periods,inf_imp(:,i+2,1)','--b','linewidth',2)
    plot(1:n_periods,inf_imp(:,i+3,1)',':b','linewidth',2)
    plot(1:n_periods,inf_imp(:,i+4,1)','-.b','linewidth',2)
    title('\textbf{Inflation}','Interpreter','latex')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
    % 341 goods, output
    subplot(3,1,3)
    plot(1:n_periods,mc_imp(:,i,1)','-b','linewidth',2)
    hold on
    plot(1:n_periods,mc_imp(:,i+1,1)','-b+','linewidth',2)
    plot(1:n_periods,mc_imp(:,i+2,1)','--b','linewidth',2)
    plot(1:n_periods,mc_imp(:,i+3,1)',':b','linewidth',2)
    plot(1:n_periods,mc_imp(:,i+4,1)','-.b','linewidth',2)
    title('\textbf{Real MC}','Interpreter','latex')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
    legend1=legend('case1','case2','case3','case4','case5','Location','best')
    fig341 = gcf
    fig341.PaperUnits = 'normalized';
    fig341.PaperPosition = [.1 0 .8 1];
    fig341.PaperPositionMode = 'manual';
    filename = ['IRFs_mu_parameter_params_' num2str(k) '_sectors_' num2str(n_sectors) '.eps'];
    saveas(fig341,filename); 
    filename = ['IRFs_mu_parameter_params_' num2str(k) '_sectors_' num2str(n_sectors) '.pdf'];
    saveas(fig341,filename); 